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The phase behavior of the lattice restricted primitive model (RPM) for ionic systems with addi- 
tional short-range nearest neighbor (nn) repulsive interactions has been studied by grand canonical 
Monte Carlo simulations. We obtain a rich phase behavior as the nn strength is varied. In particular, 
the phase diagram is very similar to the continuum RPM model for high nn strength. Specifically, 
we have found both gas- liquid phase separation, with associated Ising critical point, and first-order 
liquid-solid transition. We discuss how the line of continuous order-disorder transitions present for 
the low nn strength changes into the continuum-space behavior as one increases the nn strength 
and compare our findings with recent theoretical results by Ciach and Stell [Phys. Rev. Lett. 91, 
060601 (2003)]. 
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I. INTRODUCTION 

In spite of the considerable progress made in the last 
decade in getting a clear picture of critical behavior and 
phase separation in systems dominated by Coulombic in- 
teractions, ionic systems remain the subject of intense 
interest. For a symmetric 1:1 electrolyte system, for 
example, recent experiments Q, |^ H, ^ |^ and simula- 
tions ^6, 7, 8, 9, 10] strongly support three-dimensional 
Ising-like criticality as the asymptotic behavior. One of 
the most basic and successful models of ionic fluids is the 
restricted primitive model (RPM), in which the ions are 
viewed as equisized hard spheres carrying positive and 
negative charges of the same magnitude. Used for both 
theoretical and Monte Carlo simulations, the RPM model 
is able to characterize properly the vapor-liquid phase 
transition observed in electrolyte solutions [l^ Ha. n3t , 
as well the solid- liquid transitions of molten salts [Ij] . 

In recent years, in order to understand better critical- 
ity in the RPM modeLa lattice version of this model has 
been introduced [Illiilllliilliiillillllll. In 
this model, the positions of the positive and negative ions 
are restricted to the sites of an underlying lattice with a 
spacing equal to the ionic diameter. The most striking 
feature of this model is the presence of an order-disorder 
transition, which is absent in the continuous version of 
the RPM. There is no gas-liquid transition and the co- 
existence is between a low-density disordered phase and 
an antiferromagnetically ordered high-density phase; the 
transition is continuous (Neel-type line) above and first- 
order below a tricritical point. By contrast, non-ionic 
fluids have the same critical behavior as the lattice Ising 
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model such that continuum and lattice models are essen- 
tially equivalent f24| . 

Although the presence of an underlying lattice natu- 
rally favors the appearance of charge ordering, it is not 
completely obvious why the lattice and continuum mod- 
els of RPM present such a different critical behavior. A 
possible exp lanation has been advanced by Ciach and 
Stell flj, Hsl flfll ] using a formalism based on the Landau- 
Ginzburg- Wilson approach. They proposed that in con- 
trast to the uncharged systems with short-range interac- 
tions, where the long-wavelength fluctuations dominate 
and the lattice structure is irrelevant, in ionic systems 
the short-wavelength charge fluctuations are the most 
important. In this case, the short-distance properties 
of the system, such as the lattice structure or the shape 
of the short-range potentials added to the RPM model, 
become important and different phase diagrams can be 
obtained, i.e. there is no universality. In fact, Ciach and 
Stell have predicted that for a model system with addi- 
tional short-range interactions added to the RPM model 
both gas-liquid and tricritical points can be thermody- 
namically stable. Also, more recently [iM they pro- 
posed that when repulsive nearest neighbor interactions 
are included to the lattice RPM model, the phase dia- 
gram obtained should be qualitatively the same as in the 
continuum RPM model. These short-range interactions 
could represent the interaction between the ions and the 
particles of the solvent in which the ions are dissolved. 

In this paper we extend the lattice RPM model for 
ionic systems introduced in Ref. 22], where a short- 
range attractive potential was supplemented to the lat- 
tice RPM. Now, in addition to the hard-core and elec- 
trostatic interactions, some short-range repulsive inter- 
actions between the ions are included. We use grand 
canonical Monte Carlo simulations, combined with his- 
togram reweighting psf and mixed-field finite-size seal- 
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FIG. 1: Two-dimensional projection of the three-dimensional 
lattice structure used in our simulations. In addition to the 
electrostatic potential, a charged particle in site will be af- 
fected by a repulsion with the other particles on the first near- 
est neighbor sites (1, 2, 3 and 4, plus two more off the plane). 
For the second nearest neighbor sites (5, 6, 7 and 8, and the 
corresponding off-plane positions) just the electrostatic po- 
tential is considered. 

ing techniques, to obtain the coexistence curves and 
the associated critical points. The paper is organized as 
follows. The model and the computational details are 
given in Sec. ^ Results are discussed in Sec. IIIII We 
close in Sec. llVl with summary and conclusions. 

II. MODEL SYSTEM AND SIMULATION 
METHODS 

The model used here is essentially the same as that 
of Ref. , but with repulsive rather than attractive in- 
teractions. We consider a system of 2N charged hard 
spheres of equal diameter cr, half of them carrying charge 
+q and half charge —q, interacting through the pair po- 
tential 

{Ml. ifr >a 

U,j = \ ^"^v (1) 
[ -l-oo if rij < a , 

where D is the dielectric constant of the structureless 
solvent in which the ions are immersed. In addition to the 
electrostatic potential and hard-core given in Eq. we 
include a short-range repulsive potential of strength J > 
between the first (nn) nearest neighbor ions, regardless 
of their charge. For J = the lattice RPM model is 
recovered. These ions are restricted to the sites of a three- 
dimensional simple cubic lattice with a unit cell length 
I. In Fig. ^ we show a two-dimensional projection of 
the lattice structure used in our simulations. Reduced 
quantities are defined as follow 

where a is the ion diameter, V is the volume of the sys- 
tem and Eq = q'^ /Da is the Coulomb energy between two 
ions at close contact. The reduced chemical potential, /i*, 
is defined so that at the limit of high temperatures and 
low densities, fi* 2T* InNa^ /V, where the factor 2 



comes from the presence of two ions per minimal neutral 
"molecule" inserted or deleted in the simulations. Using 
these definitions, the effect of the repulsive interactions 
on the properties of the lattice RPM model can be mon- 
itored by changing the reduced energy parameter, J*. 

The simulations were performed using the discretiza- 
tion methodology introduced by Panagiotopoulos and 
Kumar ,!&] . In this approach the allowed positions for 
the centers of the ions are on a simple cubic grid of char- 
acteristic length /. Also, the lattice discretization pa- 
rameter is defined as C — ^/^ such that the lattice and 
continuum limits can be reproduced by changing C. In 
fact for C = 10 the results were nearly indistinguishable 
from the continuum model and the critical parameters 
of the RPM were well reproduced 0, 1131 . In this work 
we study the lattice RPM model corresponding to C = 1- 
The Ewald sums were performed with conducting bound- 
ary conditions, using 518 Fourier-space wave vectors and 
real-space damping parameter k = 5. 

We used grand canonical Monte Carlo (GCMC) sim- 
ulations with pair additions and removals at each time 
step. To enhance acceptance of the insertion and removal 
steps we used distance-biased sampling, introduced in 
Ref. m. Multihistogram reweighting j2l, H 113 tech- 
niques were used to analyze the simulation data. For 
the critical region we used mixed-field finite size scal- 
ing (FSS) analysis proposed by Bruce and Wilding [2^ . 
which accounts for the lack of symmetry between coexist- 
ing phases in fluids. We did not attempt to incorporate 
corrections for pressure mixing in the scaling fields, as 
any such effects are expected to be small 31J. Typical 
runs involve 2— 5 x 10^ Monte Carlo steps (MCs) for equi- 
libration and 2 — 9 X 10* MCs for production. Statistical 
uncertainties for the critical parameters were produced 
from 8 to 24 independent runs, depending on system size, 
at near critical conditions, with different seeds used for 
the random number generator routine "ran2" of Ref. [s^ ■ 



III. RESULTS AND DISCUSSION 

In this section we present our results obtained for the 
lattice RPM model with nn repulsion. Since J* — rep- 
resents the pure lattice RPM model, for which the phase 
diagram is known, we begin our discussion analyzing the 
dependence of the density on the chemical potential when 
we increase the nn repulsive strength J* . 

Figure |21 shows some isotherms calculated for differ- 
ent repulsive strength J*. For J* = 0.01, for instance, 
Fig. 121 (a) shows that the dependence of the density on 
the chemical potential is smooth for T* = 0.14, but for 
T* = 0.11 there is a discontinuity at /i* = -1.679. The 
same behavior persists until J* — 0.06. A closer in- 
spection of the configurations generated at this chemical 
potential, Fig. reveals an order-disorder phase transi- 
tion. In this regime the electrostatic interaction drives 
the phase separation and only a tricritical point is ob- 
served. The effect of the short-range repulsion can be 






FIG. 3: (color online) Coexisting phases generated for J* = 
0.01. The temperature is T* — 0.11 and the densities are (a) 
p* = 0.19 and (b) p* = 0.93. 




noticed only as a decrease of the tricritical temperature 
and increase of the corresponding density, as shown in 
Fig. ^ our estimate for the phase diagram as a function 
of the repulsive strength J*, obtained from histogram 
reweighting. We give an estimate of the location of the 
tricritical point, as shown in Fig.^ based on a linear ex- 
trapolation of the coexisting lines, as expected for d — 3 
tricriticality. We did not find any evidence of a gas-liquid 
phase separation for J* between and 0.06. 
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FIG. 2: Isotherms calculated in grand canonical simulations 
(L* = 12) for different repulsive strength, (a) J* = 0.01, (b) 
0.03, (c) 0.05 and (d) 0.06. The dashed line marks the approx- 
imated location of the order-disorder phase transition. The 
solid lines are just guides to the eye. Statistical uncertainties 
are smaller than the symbol size. 



0.125- 



* 0.1 




0.075 



0.2 



FIG. 4: Phase diagram as a function of the repulsive strength 
J*. Points (open circles) from top to bottom are for J* = 
0.01, 0.03, 0.05, 0.06, and 0.3, respectively. Th e C = 1 lattice 
RPM results of Panagiotopoulos and Kumar |lql (J* = 0), 
are shown as filled circles. Tricritical points ( x ) were obtained 
using a linear extrapolation (dotted lines) of coexistence lines, 
as expected for d = 3 tricriticality. The critical point (-I-) for 
J* = 0.3 was estimated from finite size scaling using L* — 15. 
Statistical uncertainties are smaller than the symbol size. 



When we increase the repulsive strength J* the nearest 
neighbor occupancy becomes less favorable. Consider the 
reduced electrostatic energy defined as follows: 



7-7-* _ ^ij _ ^i^j 
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FIG. 5: (color online) Configurations generated in the inter- 
mediate density region for (a) J* = 0.07 and (b) 0.1. In (a) 
the density is p* — 0.17 and the temperature is T* — 0.04016, 
while in (b) p* = 0.5 and T* = 0.06. 



where r. 
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j I a is the reduced separation between two 
ions of valences Zi and Zj . In the lattice model depicted in 
Fig.fflthe distance between the nearest neighbor sites is 
while for the second nearest neighbor ones is \f^l, where I 
is the characteristic length of the simple cubic grid we are 
using in the simulations. Since the lattice RPM model 
corresponds to / = cr, if we consider a = 1, the energy 
between two ions in the first nearest neighbor sites is 
ZiZj + J*, while for the second nearest neighbor ones is 
ZiZj/yj2. Therefore, for J* = 1 - l/^/2 w 0.3, if the site 
in Fig. n is populated by a negative ion [zi = —1), a 
positive ion [zj = +1) will be affected by the same energy 
-1/V2 and can occupy first and second nn sites with 
equal probabilities. Therefore, we expect competition 
between first and second nearest neighbor occupancy in 
the region around ,/* = 0.3 . In the limit of J* — > cx) the 
occupation of the first nn sites is prohibited. 

In fact, we have found that above J* = 0.065 the phase 
diagram starts to change. While in Fig. [21 we have only 
a first-order phase transition between a low density dis- 
ordered state and a high density antiferromagnetically 
ordered state, the simulations suggest the existence of 
structured configurations in the intermediate density re- 
gion, between the disordered and antiferromagnetically 
ordered states. In Fig. we show some of these con- 
figurations. We can see in Fig. |H| that there is a range 
of chemical potentials over which a plateau at p* = 0.5 
starts to evolve, where the configurations are similar to 
that of Fig.|5l(b). 

A run initiated either at low density disordered state 
or at a high density ordered state [see e.g. Figs.|31(a) and 
(b)] it would convert to this plateau. Although hysteresis 
loops are expected whenever first-order phase transitions 
are present, our GCMC method is unable to characterize 
these configurations properly. Configurations similar to 
Fig. [51 (b) remain apparently stable even after 1.8x10^ 
Monte Carlo steps, with a very low acceptance for the 
Monte Carlo moves (< 0.01%), due to the low temper- 
ature and high density. Such structures could represent 
metastable states, but we cannot determine this conclu- 
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FIG. 6: (color online) Isotherms calculated in grand canonical 
simulations {L* — 12) for different repulsive strength, (a) 
J* = 0.065, (b) 0.07 and (c) 0.1. Statistical uncertainties are 
smaller than the symbol size. 



sively based on our simulations. 

For J* = 0.3, on the other hand, the simulations have 
produced a much clearer picture, as shown in Fig. [3 
Once again, hysteresis loops were observed but now for 
both low and high-density regions. Also, the plateau ob- 
served at p* = 0.5 disappears completely. For the low- 
density region Fig. [71(b) shows that the phase separation 
is between two disordered phases which we identify as a 
gas-liquid phase transition. 
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FIG. 7; (color online) Density versus chemical potential for 
J* = 0.3. Statistical uncertainties are smaller than the sym- 
bol size. 



In order to characterize the critical point associated 
with the gas-liquid transition, we have used FSS |26| 
analysis, which accounts for the lack of symmetry be- 
tween coexisting phases in fluids. Briefly stated, for 
one-component systems we define an ordering parame- 
ter M = N — sU, where s is the field-mixed parameter, 
such that at criticality the normalized probability dis- 
tribution at a given system size, has a universal 
form for every fluid in a given universality class, with 
X = A{M — Mc). As stated in Section ^ we did not at- 
tempt to incorporate pressure mixing in the scaling fields. 

In Fig. IHlwe show the collapse of the measured Pl{x) 
on the universal Ising ordering operator distribution for 
the system sizes L* = 12, 15 and 18. The quality of 
the collapse on the universal three-dimensional Ising crit- 
ical distribution is better for the intermediate system 
size. We attribute the discrepancies for the larger sys- 
tem size to inadequate sampling, especially around the 
less probably states of a; = 0. For the smaller sys- 
tem size, there are too few particles in the simulation 
for adequate mapping on the universal distribution, as 
also observed previously for the RPM 0- We suggest 
that the system presents critical behavior compatible 



FIG. 8: (color online) Ordering operator distribution for = 
1 and J* = 0.3, for L* = 12 (triangles), L* = 15 (squares) 
and L* = 18 (circles). The L* = 15 and 18 curves have been 
displaced vertically for visual clarity. Lines are for the three- 
dimensional Ising universality class (data courtesy of N. B. 
Wilding). 



with Ising-like behavior, with a critical point located 
at T* ^ 0.04331 ± 0.00005, p* = 0.049 ± 0.006 and 
/i* = -1.1537 ± 0.0005 for a system size L* = 15. The 
critical temperature is reduced relative to the continuum 
RPM estimate 0, T* = 0.0489 ± 0.0003, mainly due to 
the addition of nearest neighbor repulsion in our model. 
The critical density, on the other hand, is considerable 
lower than = 0.076 ± 0.003, the continuum RPM es- 
timate 0. In Fig.El we show the gas-liquid coexistence 
curves for the system sizes used in Fig. |H1 along with 
the continuum RPM result obtained from the finely dis- 
cretized simulations of Panagiotopoulos and Kumar . 

For the high-density region Fig. [7| (a) suggests a first- 
order phase transition between a disordered high-density 
phase and an ordered state. Since these hysteresis loops 
still persist for a high temperature, we propose that the 
order-disorder phase transition with a tricritical point, 
observed for the low J* region, is replaced by a usual 
continuum-space behavior (first-order phase transition) 
as one increases J* . In Fig. ^| we show a typical con- 
figuration for this high density region. Since we are in- 
creasing the first nn repulsion, the system evolves to a 
less denser configuration, instead of the antiferromagnet- 
ically ordered state observed for the low J* domain. This 
behavior is essentially the same of the theoretical predic- 
tions of Ciach and Stell's approach (see Fig. 1 (d) of 
Ref. ^^). We did not include the high density region for 
J* = 0.3 in Fig. 21 since GCMC simulations are unable 
to produce adequately sampled states at this very high 
densities. 

From the phase diagram presented in Fig. 0] we can 
predict how the line of continuous order-disorder tran- 
sitions present for the low nn repulsion changes into 
the continuum-space behavior as one increases J* . For 
the low J* region, the gas-liquid transition remains 
metastable into the two-phase region associated with the 
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Ginzburg- Wilson field-theoretical approach they have 
found a phase diagram very similar to Fig. ^ in the limit 
of strong nn repulsion. Also, they have obtained the 
evolution of the phase diagram from order-disorder to 
a fluctuation-induced first-order charge-ordered-charge- 
disordered transition for high densities Although 
GCMC is unable to properly characterize high-density 
phases, our simulation results predict that the sequence 
of phase diagrams observed in the Ciach and Stell's model 
is (a)^(b)^(d) in Fig. 1 of Ref. when we increase 
the nn strength. 
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FIG. 9: (color online) Gas-liquid coexistence curves for the 
lattice RPM model supplemented by a nearest neighbor re- 
pulsive strength J* — 0.3, for L* — 12 (squares), L* — 15 
(triangles) and L* — 18 (diamonds). The continuum RPM 
curves (circles) were taken from the finely discretized sim- 
ulations from Panagiotopoulos and Kumar [l^. Statistical 
uncertainties smaller than the symbol size have been omit- 
ted. 




FIG. 10: (color online) Typical configuration observed in 
the high-density region for J* = 0.3. The temperature is 
T* = 0.07 and the density is p* = 0.75. 



order-disorder phase separation. When the nn strength 
is increased, the gas-liquid critical temperature increases, 
while the stable tricritical point decreases. Eventually, 
these two temperatures become of the same order, and 
the metastable critical point approaches the coexistence 
line of the order-disorder phase separation. This is essen- 
tially the behavior observed in Fig. El From this point, 
the increase of J* makes the tricritical point metastable. 

Recently Ciach and Stell [3 have proposed a model 
where the same nn repulsion was added to the lat- 
tice RPM model. Using a model based on a Landau- 



IV. CONCLUSIONS 

In summary, we have used grand canonical Monte 
Carlo simulation and histogram reweighting techniques 
to study phase transitions in a lattice RPM model where, 
in addition to the Coulomb and hard-core interactions, 
some short-range repulsive interactions between the ions 
are added to the model. Phase diagrams for different 
short-range strength have been obtained. Our simula- 
tion results reveal a phase diagram strongly dependent 
on the nn parameter J*. Specifically, for weak nn repul- 
sion and C = 1 only order-disorder phase coexistence and 
a tricritical point are observed, since the phase separa- 
tion is driven by the electrostatic interactions. Increasing 
J* makes the nearest neighbor occupancy becomes less 
favorable and the phase diagram starts to change. The 
nn exclusion on the simple cubic lattice used in our simu- 
lations prevents the order-disorder phase coexistence and 
tricriticality and for the low density region a gas-liquid 
phase separation appears. The critical point was esti- 
mated to belongs to the Ising universality class, with crit- 
ical parameters at the same order of the continuum-space 
estimates. Thus our simulation results confirm qualita- 
tively most of the theoretical predictions of Ciach and 
SteU's approach [THH. 
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